
options(scipen=20)

library(MASS)
library(WhatIf)
library(RColorBrewer)
library(simcf)
library(verification)
library(foreign)
library(nnet)
library(xtable)
library(apsrtable)
library(ggplot2)
library(dplyr)
library(stargazer)
col <- brewer.pal(3,"Dark2")

anes<-read.dta("anes.dta")
anes<- anes %>% mutate(ai=protest+petition+don_rel+don_soc+community+meet+volunteer+official)

ncps<-read.dta("ncps.dta")
ncps<-ncps %>% mutate(ai=petition+community+meeting+org+letter+donate+dem+elect)

# FIRST, STRAIGHT MODELS, and by WHITE/NONWHITE #

############
### ncps ###
############

m1<-lm(cpi~prox2+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
         incmis+polint+church+republican+independent+mode+nonwhite, 
       data=ncps,weights=weight)

wncps<-ncps[ which(ncps$nonwhite==0), ]

m2<-lm(cpi~prox2+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
         incmis+polint+church+republican+independent+mode, 
       data=wncps)

nwncps<-ncps[ which(ncps$nonwhite==1), ]

m3<-lm(cpi~prox2+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
         incmis+polint+church+republican+independent+mode, 
       data=nwncps)

##########
## ANES ##
##########

m4<-lm(disc~prox2+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
         incmis+polint+church+republican+independent+nonwhite, 
       data=anes,weights=weight)

wanes<-anes [ which(anes$nonwhite==0), ]

m5<-lm(disc~prox2+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
         incmis+polint+church+republican+independent, 
       data=wanes,weights=weight)

nwanes<-anes[ which(anes$nonwhite==1), ]

m6<-lm(disc~prox2+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
         incmis+polint+church+republican+independent, 
       data=nwanes,weights=weight)



stargazer(m1,m2,m3,m4,m5,m6, no.space=TRUE,
          single.row = FALSE,
          font.size = "footnotesize",column.sep.width="1.5pt",
          covariate.labels=c("Proximal",
                             "Personal",
                             "Efficacy",
                             "Female","Age: 18-34","Age: 65+",
                             "Education","20k-40k","40k-60k","60k-80k","80k-100k","100k+",
                             "Missing Inc","Political Int","Religiosity",
                             "Republican","Independent","Mode","Nonwhite",
                             "Constant"),
          dep.var.labels.include = FALSE,
          model.numbers=FALSE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          column.labels = c("Total","White", "Nonwhite", "Total","White","Nonwhite"),
          dep.var.labels   = "")


#### STEP 2: DISCRIMINATION and CONTACT --> PARTICIPATION #####

##########
## NCPS ##
##########

m1<-lm(ai~prox2+pers+cpi+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
         incmis+polint+church+republican+independent+mode+nonwhite, 
       data=ncps,weights=weight)

wncps<-ncps[ which(ncps$nonwhite==0), ]

m2<-lm(ai~prox2+pers+cpi+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
         incmis+polint+church+republican+independent+mode, 
       data=wncps)

nwncps<-ncps[ which(ncps$nonwhite==1), ]

m3<-lm(ai~prox2+pers+cpi+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
         incmis+polint+church+republican+independent+mode, 
       data=nwncps)

##########
## ANES ##
##########

m4<-lm(ai~prox2+pers+disc+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
         incmis+polint+church+republican+independent+nonwhite, 
       data=anes,weights=weight)

wanes<-anes [ which(anes$nonwhite==0), ]

m5<-lm(ai~prox2+pers+disc+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
         incmis+polint+church+republican+independent, 
       data=wanes,weights=weight)

nwanes<-anes[ which(anes$nonwhite==1), ]

m6<-lm(ai~prox2+pers+disc+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
         incmis+polint+church+republican+independent, 
       data=nwanes,weights=weight)

stargazer(m1,m2,m3,m4,m5,m6, no.space=TRUE,
          single.row = FALSE,
          font.size = "footnotesize",column.sep.width="1.5pt",
          covariate.labels=c("Proximal",
                             "Personal",
                             "Injustice",
                             "Discrimination",
                             "Efficacy",
                             "Female","Age: 18-34","Age: 65+",
                             "Education","20k-40k","40k-60k","60k-80k","80k-100k","100k+",
                             "Missing Inc","Political Int","Religiosity",
                             "Republican","Independent","Mode","Nonwhite",
                             "Constant"),
          dep.var.labels.include = FALSE,
          model.numbers=FALSE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          column.labels = c("Total","White", "Nonwhite", "Total","White","Nonwhite"),
          dep.var.labels   = "")


##############################################################################
#### STEP 3: Injustice/Discrimination*Proximal Contact --> PARTICIPATION #####
##############################################################################


m1<-lm(ai~prox2*cpi+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
         incmis+polint+church+republican+independent+mode+nonwhite, 
       data=ncps,weights=weight)

wncps<-ncps[ which(ncps$nonwhite==0), ]

m2<-lm(ai~prox2*cpi+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
         incmis+polint+church+republican+independent+mode, 
       data=wncps)

nwncps<-ncps[ which(ncps$nonwhite==1), ]

m3<-lm(ai~prox2*cpi+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
         incmis+polint+church+republican+independent+mode, 
       data=nwncps)

##########
## ANES ##
##########

m4<-lm(ai~prox2*disc+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
         incmis+polint+church+republican+independent+nonwhite, 
       data=anes,weights=weight)

wanes<-anes [ which(anes$nonwhite==0), ]

m5<-lm(ai~prox2*disc+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
         incmis+polint+church+republican+independent, 
       data=wanes,weights=weight)

nwanes<-anes[ which(anes$nonwhite==1), ]

m6<-lm(ai~prox2*disc+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
         incmis+polint+church+republican+independent, 
       data=nwanes,weights=weight)
summary(m6)
stargazer(m1,m2,m3, m4,m5,m6, no.space=TRUE,
          single.row = FALSE,
          font.size = "footnotesize",column.sep.width="1.5pt",
          order=c(1,2,3,22,23,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,24),
          covariate.labels=c("Proximal",
                             "Injustice",
                             "Discrimination",
                             "Prox*Injustice",
                             "Prox*Discrimination",
                             "Personal",
                             "Efficacy",
                             "Female","Age: 18-34","Age: 65+",
                             "Education","20k-40k","40k-60k","60k-80k","80k-100k","100k+",
                             "Missing Inc","Political Int","Religiosity",
                             "Republican","Independent","Mode","Nonwhite",
                             
                             "Constant"),
          dep.var.labels.include = FALSE,
          model.numbers=FALSE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          column.labels = c("Total","White", "Nonwhite", "Total","White","Nonwhite"),
          dep.var.labels   = "")


##########################
#### Make Figure 1-2 #####
##########################

model_ncps1<-ai~prox2*cpi+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
  incmis+polint+church+republican+independent+mode+nonwhite

model_ncps2<-ai~prox2*cpi+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
  incmis+polint+church+republican+independent+mode

sims <- 10000
pe<-c(m1$coefficients)
vc<-vcov(m1)
simbetas <- mvrnorm(sims,pe,vc)

xhyp <- seq(0,9,1)  
nscen <- length(xhyp)  

scen1 <- scen2 <- cfMake(model_ncps1, ncps, nscen)  
for (i in 1:nscen) {
  scen1 <- cfChange(scen1, "cpi", x = xhyp[i], scen = i)
  scen1 <- cfChange(scen1, "prox2", x = 0, scen = i)
  
  scen2 <- cfChange(scen2, "cpi", x = xhyp[i], scen = i)
  scen2 <- cfChange(scen2, "prox2", x = 1, scen = i)
}
y1 <- linearsimev(scen1, simbetas) 
y2 <- linearsimev(scen2, simbetas) 

pe<-c(m2$coefficients)
vc<-vcov(m2)
simbetas <- mvrnorm(sims,pe,vc)

scen1 <- scen2 <- cfMake(model_ncps2, wncps, nscen)  
for (i in 1:nscen) {
  scen1 <- cfChange(scen1, "cpi", x = xhyp[i], scen = i)
  scen1 <- cfChange(scen1, "prox2", x = 0, scen = i)
  
  scen2 <- cfChange(scen2, "cpi", x = xhyp[i], scen = i)
  scen2 <- cfChange(scen2, "prox2", x = 1, scen = i)
}
y3 <- linearsimev(scen1, simbetas) 
y4 <- linearsimev(scen2, simbetas) 

pe<-c(m3$coefficients)
vc<-vcov(m3)
simbetas <- mvrnorm(sims,pe,vc)

scen1 <- scen2 <- cfMake(model_ncps2, nwncps, nscen)  
for (i in 1:nscen) {
  scen1 <- cfChange(scen1, "cpi", x = xhyp[i], scen = i)
  scen1 <- cfChange(scen1, "prox2", x = 0, scen = i)
  
  scen2 <- cfChange(scen2, "cpi", x = xhyp[i], scen = i)
  scen2 <- cfChange(scen2, "prox2", x = 1, scen = i)
}
y5 <- linearsimev(scen1, simbetas) 
y6 <- linearsimev(scen2, simbetas) 


pr1<-cbind(y1$pe,y1$lower,y1$upper)
pr1<-data.frame(pr1)
pr1<-pr1 %>% mutate(contact="No Contact",
                    race="0")

pr2<-cbind(y2$pe,y2$lower,y2$upper)
pr2<-data.frame(pr2)
pr2<-pr2 %>% mutate(contact="Proximal Contact",
                    race="0")

pr3<-cbind(y3$pe,y3$lower,y3$upper)
pr3<-data.frame(pr3)
pr3<-pr3 %>% mutate(contact="No Contact",
                    race="1")

pr4<-cbind(y4$pe,y4$lower,y4$upper)
pr4<-data.frame(pr4)
pr4<-pr4 %>% mutate(contact="Proximal Contact",
                    race="1")

pr5<-cbind(y5$pe,y5$lower,y5$upper)
pr5<-data.frame(pr5)
pr5<-pr5 %>% mutate(contact="No Contact",
                    race="2")

pr6<-cbind(y6$pe,y6$lower,y6$upper)
pr6<-data.frame(pr6)
pr6<-pr6 %>% mutate(contact="Proximal Contact",
                    race="2")


d<-rbind(pr1,pr2,pr3,pr4,pr5,pr6)
colnames(d)<-c("pe", "lower", "upper","contact","race")

lab<-c("0"="Total Sample","1"="Whites","2"="Nonwhites")

x<-xhyp
df<-data.frame(x,d)
p<-ggplot(df, aes(x=x, y=pe, group=contact))
p1<-p+ylim(1.5,5.2)+geom_line(aes(linetype=contact))+geom_ribbon(data=df,aes(ymin=lower, ymax=upper,fill=contact),alpha=.25)+
  #scale_fill_manual(values=c("gray0", "gray70"))+
  facet_grid( . ~ race, labeller=as_labeller(lab))+
  theme_bw(base_size=14,base_family="Times")+
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="bottom",
        panel.grid.minor=element_line(colour=NA))+
  theme(plot.margin=unit(c(.5,1.5,.5,.5),"cm"))+
  theme(strip.background = element_rect(fill="white"))+
  labs(y="\nExpected Value of Participation\n",x="\nA sense of Injustice\n",
       title="")

ggsave(filename="NCPS_fig2.pdf", height = 4,width = 6, plot=p1,dpi=300)

#############
### ANES ####
#############

model_anes1<-ai~prox2*disc+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
  incmis+polint+church+republican+independent+nonwhite

model_anes2<-ai~prox2*disc+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+inc101+
  incmis+polint+church+republican+independent

sims <- 10000
pe<-c(m4$coefficients)
vc<-vcov(m4)
simbetas <- mvrnorm(sims,pe,vc)

xhyp <- seq(0,3,1)  
nscen <- length(xhyp)  

scen1 <- scen2 <- cfMake(model_anes1, anes, nscen)  
for (i in 1:nscen) {
  scen1 <- cfChange(scen1, "disc", x = xhyp[i], scen = i)
  scen1 <- cfChange(scen1, "prox2", x = 0, scen = i)
  
  scen2 <- cfChange(scen2, "disc", x = xhyp[i], scen = i)
  scen2 <- cfChange(scen2, "prox2", x = 1, scen = i)
}

y1 <- linearsimev(scen1, simbetas) 
y2 <- linearsimev(scen2, simbetas) 


pe<-c(m5$coefficients)
vc<-vcov(m5)
simbetas <- mvrnorm(sims,pe,vc)

scen1 <- scen2 <- cfMake(model_anes2, wanes, nscen)  
for (i in 1:nscen) {
  scen1 <- cfChange(scen1, "disc", x = xhyp[i], scen = i)
  scen1 <- cfChange(scen1, "prox2", x = 0, scen = i)
  
  scen2 <- cfChange(scen2, "disc", x = xhyp[i], scen = i)
  scen2 <- cfChange(scen2, "prox2", x = 1, scen = i)
}

y3 <- linearsimev(scen1, simbetas) 
y4 <- linearsimev(scen2, simbetas) 

pe<-c(m6$coefficients)
vc<-vcov(m6)
simbetas <- mvrnorm(sims,pe,vc)

scen1 <- scen2 <- cfMake(model_anes2, nwanes, nscen)  
for (i in 1:nscen) {
  scen1 <- cfChange(scen1, "disc", x = xhyp[i], scen = i)
  scen1 <- cfChange(scen1, "prox2", x = 0, scen = i)
  
  scen2 <- cfChange(scen2, "disc", x = xhyp[i], scen = i)
  scen2 <- cfChange(scen2, "prox2", x = 1, scen = i)
}

y5 <- linearsimev(scen1, simbetas) 
y6 <- linearsimev(scen2, simbetas) 


pr1<-cbind(y1$pe,y1$lower,y1$upper)
pr1<-data.frame(pr1)
pr1<-pr1 %>% mutate(contact="No Contact",
                    race="0")

pr2<-cbind(y2$pe,y2$lower,y2$upper)
pr2<-data.frame(pr2)
pr2<-pr2 %>% mutate(contact="Proximal Contact",
                    race="0")

pr3<-cbind(y3$pe,y3$lower,y3$upper)
pr3<-data.frame(pr3)
pr3<-pr3 %>% mutate(contact="No Contact",
                    race="1")

pr4<-cbind(y4$pe,y4$lower,y4$upper)
pr4<-data.frame(pr4)
pr4<-pr4 %>% mutate(contact="Proximal Contact",
                    race="1")

pr5<-cbind(y5$pe,y5$lower,y5$upper)
pr5<-data.frame(pr5)
pr5<-pr5 %>% mutate(contact="No Contact",
                    race="2")

pr6<-cbind(y6$pe,y6$lower,y6$upper)
pr6<-data.frame(pr6)
pr6<-pr6 %>% mutate(contact="Proximal Contact",
                    race="2")


d<-rbind(pr1,pr2,pr3,pr4,pr5,pr6)
colnames(d)<-c("pe", "lower", "upper","contact","race")

lab<-c("0"="Total Sample","1"="Whites","2"="Nonwhites")

x<-xhyp
df<-data.frame(x,d)
p<-ggplot(df, aes(x=x, y=pe, group=contact))
p1<-p+ylim(1.5,3.6)+geom_line(aes(linetype=contact))+geom_ribbon(data=df,aes(ymin=lower, ymax=upper,fill=contact),alpha=.25)+
  #scale_fill_manual(values=c("gray0", "gray70"))+
  facet_grid( . ~ race, labeller=as_labeller(lab))+
  theme_bw(base_size=14,base_family="Times")+
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="bottom",
        panel.grid.minor=element_line(colour=NA))+
  theme(plot.margin=unit(c(.5,1.5,.5,.5),"cm"))+
  theme(strip.background = element_rect(fill="white"))+
  labs(y="\nExpected Value of Participation\n",x="\nPerceived Discrimination\n",
       title="")

ggsave(filename="ANES_fig2.pdf",height = 4,width = 6, plot=p1,dpi=300)

#####################################
########## MEDIATION ANALYSIS #######
#####################################
library(mediation)

# NCPS 
model.m<-lm(cpi~prox2+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+
              inc101+incmis+polint+church+republican+independent+mode+nonwhite, 
            data=ncps,weights=weight)
model.y<-lm(ai~prox2+pers+cpi+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+
              inc101+incmis+polint+church+republican+independent+mode+nonwhite, 
            data=ncps, weights=weight)
#out.0<-mediate(model.m, model.y, sims=1000,boot=TRUE, treat="prox2", mediator="cpi")
out.0<-mediate(model.m, model.y, sims=5000, treat="prox2", mediator="cpi")

summary(out.0)
sens.tot.ncps <- medsens(out.0, rho.by = 0.1,effect.type = "indirect", sims = 1000)
summary(sens.tot.ncps)

pdf (file = "sensplot_ncps_tot.pdf", height = 5, width = 5)
par(family = 'Times')
plot(sens.tot.ncps, sens.par = "rho", main = "Injustice",xlim=c(-.5,.5), ylim=c(-.5,.5))
dev.off()

# white
model.m<-lm(cpi~prox2+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+
              inc101+incmis+polint+church+republican+independent+mode, 
            data=wncps)
model.y<-lm(ai~prox2+pers+cpi+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+
              inc101+incmis+polint+church+republican+independent+mode, 
            data=wncps)

out.0<-mediate(model.m, model.y, sims=5000, treat="prox2", mediator="cpi")
summary(out.0)
sens.white.ncps <- medsens(out.0, rho.by = 0.1,effect.type = "indirect", sims = 1000)
summary(sens.white.ncps)

pdf (file = "sensplot_ncps_white.pdf", height = 5, width = 5)
par(family = 'Times')
plot(sens.white.ncps, sens.par = "rho", main = "Injustice",xlim=c(-.5,.5), ylim=c(-.5,.5))
dev.off()


# nonwhite
model.m<-lm(cpi~prox2+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+
              inc101+incmis+polint+church+republican+independent+mode, 
            data=nwncps)
model.y<-lm(ai~prox2+pers+cpi+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+
              inc101+incmis+polint+church+republican+independent+mode, 
            data=nwncps)

out.0<-mediate(model.m, model.y, sims=5000, treat="prox2", mediator="cpi")

summary(out.0)
sens.nw.ncps <- medsens(out.0, rho.by = 0.1,effect.type = "indirect", sims = 1000)
summary(sens.nw.ncps)



pdf (file = "sensplot_ncps_nonwhite.pdf", height = 5, width = 5)
par(family = 'Times')
plot(sens.nw.ncps, sens.par = "rho", main = "Injustice",xlim=c(-.5,.5), ylim=c(-.5,.5))
dev.off()

############################
######### ANES #############
############################

# total pop 
model.m<-lm(disc~prox2+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+
              inc101+incmis+polint+church+republican+independent+nonwhite, 
            data=anes,weights=weight)
model.y<-lm(ai~prox2+pers+disc+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+
              inc101+incmis+polint+church+republican+independent+nonwhite, 
            data=anes, weights=weight)

out.0<-mediate(model.m, model.y, sims=5000,treat="prox2", mediator="disc")

summary(out.0)
sens.tot.anes<- medsens(out.0, rho.by = 0.1,effect.type = "indirect", sims = 1000)
summary(sens.tot.anes)


pdf (file = "sensplot_anes_tot.pdf", height = 5, width = 5)
par(family = 'Times')
plot(sens.tot.anes, sens.par = "rho", main = "Perceived Discrimination",xlim=c(-.5,.5), ylim=c(-.5,.5))
dev.off()

# white
model.m<-lm(disc~prox2+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+
              inc101+incmis+polint+church+republican+independent, 
            data=wanes,weights=weight)
model.y<-lm(ai~prox2+pers+disc+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+
              inc101+incmis+polint+church+republican+independent, 
            data=wanes, weights=weight)
out.0<-mediate(model.m, model.y, sims=5000,treat="prox2", mediator="disc")

summary(out.0)
sens.white.anes<- medsens(out.0, rho.by = 0.1,effect.type = "indirect", sims = 1000)
summary(sens.white.anes)


pdf (file = "sensplot_anes_white.pdf", height = 5, width = 5)
par(family = 'Times')
plot(sens.white.anes, sens.par = "rho", main = "Percieved Discrimination",xlim=c(-.5,.5), ylim=c(-.5,.5))
dev.off()

# nonwhite
model.m<-lm(disc~prox2+pers+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+
              inc101+incmis+polint+church+republican+independent, 
            data=nwanes,weights=weight)
model.y<-lm(ai~prox2+pers+disc+inteff+female+young+old+edu4+inc40+inc60+inc80+inc100+
              inc101+incmis+polint+church+republican+independent, 
            data=nwanes, weights=weight)
out.0<-mediate(model.m, model.y, sims=5000,treat="prox2", mediator="disc")

summary(out.0)
sens.nw.anes<- medsens(out.0, rho.by = 0.1,effect.type = "indirect", sims = 1000)
summary(sens.nw.anes)


pdf (file = "sensplot_anes_nonwhite.pdf", height = 5, width = 5)
par(family = 'Times')
plot(sens.nw.anes, sens.par = "rho", main = "Perceived Discrimination", xlim=c(-.5,.5), ylim=c(-.5,.5))
dev.off()

########################
#### Make figure 3 #####
########################

# NCPS 

var<- rbind(1	,
            2	,
            3	,
            1	,
            2	,
            3	,
            1	,
            2	,
            3	)

# copied and pasted from mediation summary stats

est<-rbind(0.0648	, #total
           0.5464	,
           0.6112	,
           
           0.0248	, #white
           0.2504	,
           0.2752	,
           
           0.18379	, #nonwhite
           0.47027	,
           0.65406	)

lower<-rbind(0.0062	, #total
             0.2378	,
             0.2971 	,
             
             -0.0190	, #white
             -0.1964	,
             -0.1825	,
             
             0.05491	, #nonwhite
             -0.00239	,
             0.1995	)

upper<-rbind(0.14	, #total
             0.86	,
             0.93	,
             
             0.09	, #white
             0.71	,
             0.73	,
             
             0.34	, #nonwhite
             0.96	,
             1.13		)

group<-rbind(1	,
             1	,
             1	,
             
             2	,
             2	,
             2	,
             
             3	,
             3	,
             3)

d<-cbind(var, est,upper,lower,group)
d<-data.frame(d)
variable<-c("03","02","01")
dncps<-cbind(variable,d)
dncps<-dncps%>% mutate(data="1")

#ANES
var<- rbind(1	,
            2	,
            3	,
            1	,
            2	,
            3	,
            1	,
            2	,
            3	)

# copied and pasted from mediation summary stats

est<-rbind( 0.0375 	, #total 
            0.3875	,
            0.4250	,
            
            0.01905	, #white
            0.27043,
            0.28949,
            
            0.0653, #nonwhite
            0.6600,
            0.7253	)

lower<-rbind(0.0159	, #total 
             0.2255	,
             0.2656	,
             
             -0.00245 	, #white
             0.07779	,
             0.09921	,
             
             0.0108	, #nonwhite
             0.3383	,
             0.4079	)

upper<-rbind( 0.06	, #total
              0.55	,
              0.58	,
             
              0.04	, #white
              0.46	,
              0.48	,
             
              0.13	, #nonwhite
              0.98 	,
              1.04 		)

group<-rbind(1	,
             1	,
             1	,
             
             2	,
             2	,
             2	,
             
             3	,
             3	,
             3)

d<-cbind(var, est,upper,lower,group)
d<-data.frame(d)
variable<-c("03","02","01")
danes<-cbind(variable,d)
danes<-danes%>% mutate(data="2")

df<-rbind(dncps,danes)
colnames(df)<-c("variable","var","pe", "lower", "upper", "group","data")
df$group<-as.factor(df$group)
levels(df$group)
df$group<-factor(df$group, levels=rev(levels(df$group)))

lab<-c("1"="NCPS","2"="ANES")

p1<-ggplot(data = df, aes(x = variable, y = pe, 
                          ymin=lower, ymax = upper, group=factor(group),color=factor(group),
                          order=-as.numeric(group)))+
  scale_color_discrete(name="group",
                       h=c(0,360,180),c=75,l=40,
                       breaks=c("1","2","3"),
                       labels=c("Total","White","Nonwhite"))+
  facet_grid(.~data, labeller = as_labeller(lab))+
  geom_pointrange(size=.9,position=position_dodge(width=.5))+
  theme_bw(base_size=20,base_family="Times")+
  coord_flip()+scale_x_discrete(breaks=c("03","02","01"),
                                labels=c("Mediation Effect","Direct Effect","Total Effect"))+
  theme(panel.grid.major.x = element_blank(),panel.grid.major.y = element_blank())+
  theme(strip.background = element_rect(fill="white"))+
  theme(strip.text.x = element_text(size = 24))+
  theme(legend.title=element_blank(),
        legend.key.size=unit(2.0, 'lines'))+
  geom_hline(yintercept = 0, lty = 2)+
  theme(plot.margin=unit(c(1,2,1,1),"cm"))+
  labs(y=NULL,x=NULL,title=NULL)

ggsave(filename="NCPS_ANES_mediation.pdf", plot=p1,dpi=300,width = 30, height=20, units="cm")



